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We use uniform matrix product states (MPS) to study the (1+1)D 0(2) and 0(4) rotor models, 
which are equivalent to the Kogut-Susskind formulation of matter-free nonabelian lattice gauge 
theory on a “hawaiian earring” graph for U(l) and SU( 2), respectively. Applying tangent space 
methods to obtain ground states and determine the mass gap and the /? function, we find excellent 
agreement with known results, locating the BKT transition for 0(2) and successfully entering the 
asymptotic weak-coupling regime for 0(4). To obtain a finite local Hilbert space, we truncate in 
the space of generalized Fourier modes of the gauge group, comparing the effects of different cutoff 
values. We find that higher modes become important in the crossover and weak-coupling regimes 
of the nonabelian theory, where entanglement also suddenly increases. This could have important 
consequences for TNS studies of Yang-Mills on higher dimensional graphs. 


I. INTRODUCTION 

Nonabelian gauge theories describe the interactions re¬ 
sponsible for most of the matter we experience in our 
everyday lives. In particular, they explain the hadrons 

bound states of quarks — which include neutrons and 
protons among their most famous examples [1]. Curi¬ 
ously, the quarks inside hadrons behave as free parti¬ 
cles for the purposes of high-energy scattering (asymp¬ 
totic freedom), yet they are never observed in isolation 
(confinement [2]). These properties are also present in 
matter-free nonabelian gauge theory (pure Yang-Mills 
theory). This apparently simple theory, despite its huge 
symmetry group of local gauge transformations, resists 
exact solution and must so far be approached with ap¬ 
proximate methods such as perturbation theory and nu¬ 
merical tools such as Monte Carlo sampling [3], albeit 
with convincing successes, such as the determination of 
Hadron masses using lattice simulations [4] . 

Monte Carlo techniques are also extremely useful in 
condensed matter physics and advances have benefited 
both fields. However, in recent decades new, highly 
general techniques have arisen in condensed matter and 
quantum information that open up whole new avenues of 
numerical investigation. These techniques exploit Ten¬ 
sor Network States (TNS) [5], which efficiently repre¬ 
sent many-body states with limited entanglement. The 
best-known example is the density-matrix renormaliza¬ 
tion group (DMRG) [6] , which can be viewed [7] as a vari¬ 
ational algorithm applied to one-dimensional TNS, also 
known as Matrix Product States (MPS) [8]. These meth¬ 
ods are inherently free of the sign problem that plagues 
Monte Carlo sampling [9] and offer themselves naturally 
to simulation of real-time dynamics. 

In the last years, as TNS techniques have advanced 
(higher dimensions, more sophisticated networks, im¬ 
proved numerical tools) [10], efforts have increased to 
transfer their successes in condensed matter to quantum 
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Figure 1. Illustration of lattice gauge theory (a) on a 
“hawaiian earring” and (b) the same theory visualized dif¬ 
ferently as living on the surface of a 3D object. The Hamilto¬ 
nian is the Kogut-Susskind [29] formulation of lattice gauge 
theory for a gauge group G — example plaquette operators 
are shown in blue. The Hilbert space H (including non¬ 
physical states) is made up of systems living on the (black) 
edges Hedge = L 2 (G). For the gauge groups G = U( 1) and 
G = SU( 2) this model is equivalent to the (l+l)-dimensional 
quantum rotor model [30] for the rotation groups 0(2) and 
0(4), respectively. 


field theory, particularly with an eye toward nonabelian 
gauge theory. Important steps in this direction include 
ground state, real-time, and finite-temperature simula¬ 
tions of <p A theory [11, 12], the Schwinger model [13- 
16], SU( 2) gauge theory with matter [17] in (1+1)D, and 
quasi-one-dimensional abelian gauge theories [18], all us¬ 
ing MPS or DMRG, as well as proposals for representing 
lattice gauge theory states in higher dimensions, with a 
view toward numerics as well as analytics [19-23]. The 
tensor renormalization group (TRG) algorithm [24] has 
also been applied to q i > 4 theory [25], the Schwinger model 
[26], and the 0(2) and 0(3) models [27, 28]. 


A. This work 

Here we present an MPS study of the 0(2) and 0(4) 
quantum rotor models in 1 + ID, which are equivalent, 
respectively, to the U( 1) and SU(2) principal chiral field 
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(PCF) models. The PCF model is in turn equivalent to a 
pure lattice gauge theory on a “hawaiian earring” graph 
in the Kogut-Susskind formulation (see Figure 1), insofar 
as we do not restrict to gauge-invariant states. 

The Hamiltonians possess a global gauge-group sym¬ 
metry rather than a local gauge symmetry, but never¬ 
theless have a lot in common with Yang-Mills on more 
sophisticated graphs. Most importantly, the 0(N > 2) 
models are known to possess a single, gapped phase end¬ 
ing at the weak-coupling limit g — > 0 [30]. This is also 
observed in simulations of (3 + 1)D nonabelian lattice 
gauge theory, in which the gapped phase is confining 
[2]. In contrast, the 0(2) model has a phase transition 
(of Berezinskii-Kosterlitz-Thouless (BKT) type [31]) at 
finite coupling, transitioning into a deconfined, gapless 
phase at weaker couplings. 

The continuum limit of the rotor models, the so-called 
O(N) nonlinear sigma model [30], can be solved using 
the Bethe ansatz for N > 2 making the lattice weak- 
coupling scaling of the mass gap computable [32], The 
O(N) model has also been thoroughly investigated using 
strong-coupling expansions [33, 34], which operate on the 
same 1 + ID Hamiltonian model we study here, as well 
as high-temperature expansions (for example [35]) and 
Monte Carlo numerics ([36-39] is an incomplete selec¬ 
tion) applied to the 2D classical O(N) model. Lanczos 
diagonalization with finite-size-scaling has also been used 
[40]. In this work, we use uniform MPS to represent infi¬ 
nite, translation invariant states, applying the nonlinear 
conjugate gradient method [12] to obtain ground states 
and the MPS tangent space as an ansatz for low energy 
excitations [41, 42], determining the mass gap and the (3 
function at finite couplings and thus obtaining the phase 
diagram. 


II. THE MODEL 

A. Kogut-Susskind Hamiltonian 

The Kogut-Susskind Hamiltonian [29] on the “hawai¬ 
ian earring” graph is given by 

/— 2 + 00 r\ j— +00 

H ks(9) = ^L J 2 E 2 k -^ J 2 Re(trK4 +1 )), 

k——oo k——oo 

(i) 

where g is the coupling, a the lattice spacing, and g an 
anisotropy parameter required to ensure the renormal¬ 
ized theory is Lorentz-invariant in the continuum limit 
[43]. The Hilbert space is the tensor product of spaces 

H k = L 2 (G) 

assigned to each edge e in the graph and G is the gauge 
group. We define A a to be the Hermitian generators of 
G (for SU(2) these are the Pauli matrices A a = \a a 
with a = 1,2,3, for C7(l) there is only one A = 1). The 


operator 

El = Y. E lk 

cx. 

is the quadratic Casimir operator representing the ki¬ 
netic energy within the gauge group at edge k. The E a 
represent the infinitesimal group action 

E a . d e L e ieA a | e _Q > 

where L x implements rotations from the left, acting on 
a “position” basis as 

L x \v) = \xv) 

for x, v £ G. The tty are gauge group position operators 
defined as 

Uij\v) = t(v)ij\v), 

with t(v) an irrep of G (we choose e' e for t/(l) and the 
spin-half representation for SU(2)). This results in the 
commutator 

[E a ,Uij\ = A a,ij'Ujij. 
o' 

B. Quantum Rotor Hamiltonian 

The model (1) is known to be equivalent to a chain of 
coupled O(N) rotors (see e.g. [43]), given by 

/ +00 /— +oo 

Hk{9) = ^ 1 E E "*•"*+!. (2) 

zo, ga 

k ——oo k——oo 

where is a ./V-dimensional unit vector representing 
the fcth rotor and J k is the rotor kinetic energy. The 
normalization of J 2 is chosen to match [33]. The Hamil¬ 
tonian (2) is manifestly invariant under a global O(N) 
symmetry. The relations which underlie Hks = Hr are 
given in Table I. 


H r ,0(N) 

Hks,G = 1/(1) 

Hks, G = SU(2) 

N 

2 

4 

9 

g 2 /V 2 

g 2 / 4 

J 2 

E 2 

4 E 2 


(Re(u), Im(u)) M 

—itr(A M w) 


g = l,2 

g = 0 ... 3, Ao = 


Table I. This table shows how quantities in the rotor model (2) 
must be set to obtain equivalence to the Kogut-Susskind 
model (1) for the gauge groups G = 1/(1) and G = SU( 2). 

The continuous position basis \v) does not lend itself to 
use with MPS numerics, for which we require a discrete, 
finite basis = C d . Instead we make use of the generalized 
Fourier basis given by the Peter-Weyl theorem [44], in 
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which the matrix elements of the irreducible representa¬ 
tions (irreps) of a compact Lie group G label generalized 
Fourier modes \ij)i (the matrix element i,j of irrep l). 
The kinetic term E 2 is diagonal in this basis, with 

E 2 \ij), = i(i + i)\ij)t, 

for SU(2) and 

E 2 \n) = n 2 \n), n £ Z 

for 17(1). In the strong coupling regime g 2 1, the 
E 2 term strongly penalizes higher irreps, so we can ne¬ 
glect them to good approximation at larger g 2 , expecting 
them to become more relevant as we near weak coupling. 
Importantly, truncating the basis at a certain irrep level 
(Fourier mode) does not prevent representation of states 
invariant under the global gauge-group symmetry (since 
rotations do not mix irreps). 


III. NUMERICAL METHODS 


The uniform MPS variational class consists of states 


d 

|*(A)> = £4 

s =1 


+M 

A Sk 


_k——M 


Vr\S-M ■ ■ ■ So ■ ■ ■ S+m), 


where A s is a D x D matrix and d the dimension of the 
chosen local basis and we make the chain length infinite: 
M —> oo. For a given basis element s = (s_m, •.., s+m), 
the matrix product in square brackets determines the 
value of ( S-m ■ ■■ So ... s+mI’F^)), hence the name “ma¬ 
trix product states.” Two key features of MPS are the ef¬ 
ficiency with which quantities such as expectation values 
of local operators and correlation functions can be com¬ 
puted (requiring 0(D 3 ) multiplications), and the fact 
that the restriction to MPS form serves only to limit the 
amount of entanglement that can be present in the state. 
The dimension D is called the bond dimension and serves 
to control the degree of spatial correlations, placing an 
upper bound S < log D on the entanglement entropy, for 
example. For more background, see for example [5]. 

To make the limit M —> oo behave appropriately, we 
require A such that the transfer operator 


E := ^ A s ( 8 > A s 

S 


has spectral radius p{E) = 1 with a unique eigenvalue of 
largest magnitude (injectivity) equal to one. With this 
condition, the boundary vectors vl and vr drop out of 
all relevant calculations and (’F(A)I’I'(A)) = 1. For more 
details on using infinite, uniform MPS, see [10]. 

We set d to accommodate the dimensions of all Fourier 
modes up to a cutoff. With U( 1) all irreps are one¬ 
dimensional and we may label Fourier modes as n € Z, 
so that a cutoff is given by |n| < ?z max and 

d — 2 n max -f- 1 . 


For SU(2) we must set 

^max ^max 

d =Yl dim W) 2 = ^2( 2l + i) 2 . 

(=0 z=o 

with l = 0 , 5 , 1 , ... . 

In this study, we use values of n max up to 10 and Z max 
up to 2. The former requires d = 21, while the latter 
implies d = 55, which is unusually high for MPS numer¬ 
ics. In the algorithms of [12, 41, 45] the cost of com¬ 
putations involving nearest-neighbor operators, such as 
the potential term u e u' e+1 , scales as 0(d 4 ). We reduce 
this to 0(d 2 m), where m is the number of terms in the 
tensor product decomposition, by implementing them as 
two-site matrix product operators [7] of dimension m , 
where m < 4 for our purposes. We further accelerate 
our implementation by computing the iterands of itera¬ 
tive parts using general purpose graphics processing units 
(GPGPU’s) [46], 

With these optimizations, we apply the nonlinear con¬ 
jugate gradient (CG) method to obtain ground states 
[12], with the time-dependent variational principle [45] 
used in a pre-optimization step to provide good start¬ 
ing points for the CG algorithm. We converge all states 
up to \\Pt { a)H\^{A)}\\ < 10 -8 , where Pt{A) projects 
the energy gradient vector onto the MPS tangent space 
at 14'(A)). We then obtain low-lying excited states us¬ 
ing the methods of [41], always operating directly in the 
space of infinite, uniform MPS. All MPS algorithms used 
here are implemented as part of the open source evoMPS 
project [47]. 


IV. RESULTS 


A. Phases observed 

Since our choice of truncated basis is most appropri¬ 
ate at strong coupling, we study the system starting at 
1/g —> 0 and then approach weak-coupling as far as pos¬ 
sible, whilst maintaining accuracy. We find, for both 
the 0(2) and the 0(4) rotor, that the MPS approx¬ 
imate ground state breaks the global O(N) symmetry 
for 1/g < 1/<7sb for constant, finite D. The symmetry¬ 
breaking, for the values of D in use, is confined to a 
relatively narrow region of parameter space. Since the 
breaking of a continuous symmetry is forbidden by the 
Mermin-Wagner theorem [48], this must be a symptom 
of finite-entanglement effects [49]: The bond dimension 
needed to accurately represent the symmetric state must 
suddenly grow as we approach weak coupling. 

The 0(2) rotor is known to possess a gapless phase at 
weak coupling, characterized by algebraically decaying 
correlations, such that the correlation length is infinite 
[30]. As defined above, a uniform MPS would require 
D —► oo to accurately represent such a ground state, 
thus explaining nonphysical symmetry-breaking at finite 
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Figure 2. The location gf B of MPS ground state symmetry 
breaking (SB) as a function of the bond dimension D. We 
plot the Kogut-Susskind coupling <?| B rather than gsB to aid 
comparison. Values of g$ B were found using bisection up to 
a precision in g~ 2 of ±0.005. For 17(1) ~ 0(2), <?g B tends 
toward a finite value. This value should correspond to the 
location of the BKT transition [30]. In terms of the 0(2) rotor 
parameter the fit gives g B kt = 1.119 ± 0.004. For SU(2) ~ 
0(4), the transition does not converge for the data available. 
This is consistent with it occurring at g = g = 0. 


D. The existence of a phase transition at finite g also 
explains the narrowness of the region where symmetry¬ 
breaking begins. We expect the symmetry-breaking lo¬ 
cation gsB(D) to converge to the location of the phase 
transition as D —>• oo and indeed this convergence can be 
seen in Figure 2, where the extrapolated transition point 
g = 1.119±0.004 agrees well with a known estimate from 
strong-coupling expansions of the mass gap g ~ 1.12 [33] 
and less well with an estimate g ss 1.05 based on Pade ap¬ 
proximations of /3-functions from strong-coupling expan¬ 
sions [34] , as well as a number of other methods [50] that 
indicate g ~ 1.00 (the parameter given in these studies 
is usually x = 2/g 2 ). It is worth noting that estimating 
the transition point of a BKT transition is notoriously 
difficult due to the exponential scaling of the mass gap 
near the transition [31] and it is possible that our esti¬ 
mate would shift given data at larger bond dimensions, 
or by the use of more reliable indicators than the onset 
of nonphysical symmetry-breaking. An accurate deter¬ 
mination of <)bkt is, however, beyond the scope of this 
work. 

Despite the impossibility of representing the ground 
state precisely in the gapless phase at weak coupling, the 
scaling of von Neumann entropy and correlation length in 
MPS ground states with a range of finite D can be used 
to estimate the central charge c of the conformal field 
theory (CFT) describing the phase [51, 52], We fit data 
for D = 22, 28,34,..., 80 at \/{gyj2) = 0.75, 0.8,0.85, 0.9 
and find c = 0.992 ± 0.009, matching the known result 
of c = 1 for the 2D classical XY model [53], which is 
identical with the classical 0(2) rotor. 

We now turn to the 0(4) rotor, which is known to exist 
in a single, gapped phase down to the weak-coupling limit 
g —> 0 [30]. Here, we expect our choice of basis to become 


increasingly bad as we approach weak-coupling, due to 
the occupation of higher Fourier modes. We also expect 
greater entanglement in the exact ground state as the 
potential term, coupling nearest-neighbor edges, begins 
to dominate, and the lattice correlation length grows. 
This is not enough, however, to explain the very sudden 
occurrence of nonphysical symmetry-breaking. This is 
likely due to the “crossover” phenomenon, a property of 
the 0(N > 2) models and of nonabelian gauge theories, 
referring to persistence of strong-coupling behavior up to 
a certain region of parameter space, where weak-coupling 
behavior rapidly takes over. Despite the sudden onset of 
the weak-coupling regime, we still expect the nonphysical 
symmetry-breaking transition to disappear as D —> oo, 
as we indeed observe in Figure 2. 


B. Mass gap 

Our next source of information is the mass gap. Here, 
we can directly compare our results with the results 
of 8th-order and 6th-order strong-coupling series expan¬ 
sions (SCE) for the 0(2) and 0(4) models, respectively 
[33]. We find excellent agreement for both models up to 
the vicinity of the 0(2) phase transition and the 0(4) 
crossover region. Moving closer, Figure 3 shows that the 
mass gap descends toward zero at a finite coupling for 
0(2), whereas for 0(4) the log-linear plot shows linear 
behavior, indicating a finite mass gap for all finite cou¬ 
plings. For comparison, we plot the exact asymptotic 
weak-coupling scaling for 0(4) [32], taking into account 
speed-of-light renormalization effects due to the stark 
space-time asymmetry of the Hamiltonian discretization 
[43]. We find very good agreement with the weak- 
coupling prediction, showing that we are successfully en¬ 
tering the asymptotic scaling regime, although we also 
see from the plot that finite entanglement effects start to 
limit the accuracy (the D = 140 curve remains accurate 
for longer than the D = 91 curve for Z max = 2), as indeed 
does the Fourier mode truncation (the Z max = 2 curve is 
more accurate than the Z max = 3/2 curve for D = 140). 

In the 0(4) case it is also interesting to note that, hold¬ 
ing the Fourier cutoff at l < 2, increasing the bond dimen¬ 
sion appears to interpolate between the SCE result and 
the weak-coupling result. The effect is even clearer when 
more values of D are considered. This makes sense if we 
recall that the SCE for the mass gap perturbs the g —> oo 
ground state and first excited state, both product states 
without entanglement in the Fourier basis, by repeat¬ 
edly applying the nearest-neighbor term in the Hamilto¬ 
nian up to some order [29, 33]. The higher the order, 
the less local the correlations introduced into the states 
will be. In the same way, raising the bond dimension 
of an MPS ground state utilizing the Fourier basis allows 
longer-range correlations to be represented. In this sense, 
our MPS methods and SCE’s are very similar techniques, 
and it is not surprising that the MPS results at smaller 
bond dimensions match low-order SCE results well. 
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Figure 3. The MPS mass gap for (a) the 0(2) rotor and (b) 
the 0(4) rotor for bond dimensions D and Fourier cutoffs in 
|n| and l respectively. The strong coupling expansion of [33] 
is shown (SC), as is the weak coupling result (WC) for the 
0(4) case, which is known exactly [32]. In (b), the curves are 
adjusted by an anisotropy parameter yd) to account for the 
renormalization of the speed of light [43] (for 0(2), g is set 
to one). Near the phase transition for 0(2), and as we enter 
the weak coupling regime for 0(4), finite entanglement effects 
and, for 0(4), Fourier cutoff effects become important. 


Using similar reasoning, we can understand why the 
SCE results and the lower-Z? MPS results underestimate 
the mass gap in the 0(4) weak-coupling regime. Given 
that the ground state in this region consists of highly 
non-local structures (Wilson loop excitations of various 
sizes) [29], limiting the order of the SCE or restricting the 
amount of entanglement in the MPS should both work 
against achieving these low-energy configurations, result¬ 
ing in an overestimation of the ground state energy. In¬ 
deed, we observe significant differences on the order of 
10~ 2 in the ground state energy with D as we enter the 
crossover regime. Assuming the first excitation is repre¬ 
sented relatively accurately, this explains the underesti¬ 
mation of the gap. 

We note here that using a symmetric tensor network 
ansatz [54] might significantly extend the range of acces¬ 
sible effective bond dimensions and so enable further pen¬ 
etration into the 0(4) weak-coupling regime, although it 
would not allow access to the lowest-lying excitations of 
the 0(4) model, which break the 0(4) symmetry. For a 
model with truly local gauge symmetry, methods such as 
that of [16] are required. 


C. Beta functions 


Using the mass gap and its first derivative in the cou¬ 
pling, one can calculate (see, for example [33]) the /3 func- 


Figure 4. Beta functions determined from the MPS mass gap 
for (a) the 0(2) rotor and (b) the 0(4) rotor, together with a 
Pade approximant based on a strong-coupling expansion and 
the weak coupling result (WC) for the nonabelian case [33]. 
The estimate g bkt from Figure 2 of the location of the 0(2) 
phase transition is marked in (a) using a black rectangle. For 
0(4) it is clear from the D = 140, Z max = 3/2 curve that the 
numerical results begin to qualitatively follow the weak cou¬ 
pling behavior. However, there are clearly systematic errors 
present. This is expected because the beta function involves 
the numerical derivative of the mass gap, making it sensitive 
to small inaccuracies due to finite entanglement and Fourier 
mode truncation. 


tion as 


-mra 


(, 4 F'(g) \ - 1 

l 9 2 F(g) ) 


where F(g) = 2a(E 1 (g) - E 0 (g))/g. Using finite- 
differences to compute F ', we may compute f) functions 
from our mass gap results. 

It is also possible to use an SCE for the mass gap to 
construct a Pade approximant for the /3 function. Fur¬ 
thermore, the 0(N > 2) weak-coupling behavior of /3(g) 
is known from perturbation theory to be 


-/J(S) = (JV-2)f[ + (iV-2)jf, 


allowing this information to be incorporated, resulting 
in an approximate f3 function for all couplings for the 
0(N > 2) rotor [33]. 

We compare the Pade approximants of [33] with our 
MPS results in Figure 4, observing excellent agreement 
at stronger couplings, with the numerical results devi¬ 
ating from the approximate curve as we near the phase 
transition. In the case of 0(2), the numerical data ap¬ 
pears to predict a higher value for the phase transition 
location than the Pade approximant, in good agreement 
with our result from Figure 2. However, we also observe 
a shift in the results as the bond-dimension changes, with 
the higher-!? results corresponding to a smaller predic¬ 
tion for Obkt- This supports the possibility mentioned in 








6 


the previous section that using higher bond-dimensions 
would result in a better correspondence with the majority 
of literature results. 

The 0(4) data ceases to follow the Pade curve as we 
enter the crossover region, but does not succeed in fol¬ 
lowing the weak-coupling result accurately either. This is 
not unexpected, as both approximations are likely inac¬ 
curate in the crossover region. We do, however, see large 
variations with D and l max , particularly as we near the 
nonphysical symmetry-breaking transition. That errors 
are more visible for the /3 function than for the mass gap 
is expected since the numerical derivative amplifies small 
errors in the mass gap. We would need to reach higher 
bond dimensions and Fourier mode cutoffs to achieve ac¬ 
curate results further into the weak-coupling regime. A 
further way of reducing noise would be to compute the 
derivative F'(x) analytically from the MPS excited state. 


V. CONCLUDING REMARKS 

We have shown that tensor network state (TNS) meth¬ 
ods, in this case uniform matrix product states (MPS), 


can successfully represent states of the nonabelian quan¬ 
tum rotor model into the weak-coupling regime. The 
finite local basis, achieved through Fourier-mode trunca¬ 
tion successfully and efficiently captures strong-coupling 
physics, but becomes a more severe limitation at weak 
couplings where, additionally, the spatial entanglement 
grows substantially. 

This is promising for TNS approaches to pure non¬ 
abelian gauge theory, which is believed to possess a very 
similar phase diagram to the O(N) rotor models and, on 
the “hawaiian earring” graph, is indeed equivalent to the 
rotor models studied here. Our study also shows that 
high spatial entanglement is a feature of the theory from 
the crossover region onward, into weak-coupling. This 
may pose a challenge for numerical approaches if it car¬ 
ries over to higher dimensional nonabelian lattice gauge 
theory, since large bond dimensions may be needed to 
access the asymptotic scaling regime. 
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